% factorize A to W*H'
% varargin
% 'k'                 : rank of W & H
% 'max_iter'          : maximum number of iterations
% 'epsilon'           : additive value to prevent division by 0
% 'basis'             : specificy basis in 'W' or 'H' (for normalization)
function [W,H] = multiplicativeNMF(A,varargin)
    par.k = 5;
    par.max_iter = 100;
    par.epsilon = 1e-8;
    par.basis = 'H';
    par.verbose = 0;
    par = process_parameter(par, varargin{:});
    par.m = size(A,1);
    par.n = size(A,2);
    if par.verbose
        display('NMF parameters:')
        display(par);
    end
    
    W = rand(par.m, par.k);
    H = rand(par.n, par.k);
    for iter = 1:par.max_iter
        W = W .* (A*H) ./ (W * (H'*H) + par.epsilon);
        H = H .* (A'*W) ./ (H * (W'*W) + par.epsilon);
    end
    
    [W,H] = normalizeNMFFactor(W, H, par.basis); % normalize basis in H
end
